#include "TH1D.h"
#include "TRandom.h"
#include "TF1.h"
#include <vector>
#include "TCanvas.h"
#include "TMath.h"

using namespace std;

void flow_ep_test(){

	

	
		TRandom *myRandom = new TRandom(123);
		TH1D *hPhi = new TH1D("hPhi","hPhi of 1st event", 100, -4, 4);
		TH1D *hEP_input = new TH1D("EPinput","EP input", 100, -4, 4);
		TH1D *hEP_cal = new TH1D("EPcalculated","EP calculated",100, -4,4);
		TH1D *hEP_error = new TH1D("EPerror", "EP error", 100, -4, 4);


		for(int ievt=0; ievt < 10000; ievt++){

				TF1 *dNdPhi = new TF1("dNdphi", "1+2*[0]*cos(x-[1])+2*[2]*cos(2*(x-[3]))+2*[4]*abs(cos(3*(x-[3])))", -1*TMath::Pi(), TMath::Pi());

				double EP_input = myRandom->Uniform(-1*TMath::Pi()/2, TMath::Pi()/2);

				dNdPhi->SetParameter(0, 0); // v1 (ingnore here)
				dNdPhi->SetParameter(1, 0);
				dNdPhi->SetParameter(2, 0.1); // v2
				dNdPhi->SetParameter(3, EP_input); // v2 EP
				dNdPhi->SetParameter(4, 2); //v3

				vector <double> p_phi;
				double Qx =0;
				double Qy =0;
//				if (EP_input < 0 || EP_input > TMath::Pi()/1.5) continue;
				for( int ip=0; ip<1000; ip++){
						p_phi.push_back(dNdPhi->GetRandom());
						if (ievt==0) { hPhi->Fill(p_phi[ip]); }
						Qx += TMath::Cos( 2 * p_phi[ip] );
						Qy += TMath::Sin( 2 * p_phi[ip] );
				}

				double EP_calculated = TMath::ATan2( Qy , Qx  )/2 ;	
				
				hEP_input->Fill(EP_input);
				hEP_cal->Fill(EP_calculated);	
				hEP_error->Fill(EP_input-EP_calculated);

		}

		TCanvas *c1 = new TCanvas();
		c1->Divide(3,1);
		
		c1->cd(1); hEP_input->Draw();
		c1->cd(2); hEP_cal->Draw();
		c1->cd(3); hEP_error->Draw();
}
